## Replication for Figure 2
# In case of any questions, please email anselm.rink@gmail.com

rm(list=ls())       
library(foreign)
library(ggplot2)
library(gridExtra)
library(gtable)
setwd("")
data <- read.dta("")

data_c <- data[, c(1:2,63)]
data_c <- data_c[order(data_c$respondent_muslim, data_c$id) , ]
data_cn <- data_c[complete.cases(data_c), ]
Religion <- c(rep(1,288),rep(2,287))
data_cn$Religion <- Religion
time <- seq(1:575)
data_cn$time <- time
df <- data_cn
str(df)
df$y <- df$q27_radicalization
rm(data, data_c, data_cn, time, Religion)


#### 0.2 nbinom, simulate ####

library(MASS)
mod_nb <- glm.nb(y ~ time * Religion, data = df)
period = 288
mod_nb1 <- glm.nb(y ~ time * cos(time/period*2*pi)*sin(time/period*2*pi)*Religion, data = df)


set.seed(42)
df_pred <- stack(simulate(mod_nb1, nsim=10000))
df_pred$time <- df$time
df_pred$Religion <- df$Religion

## summarize

library(plyr)
df_pred2 <- ddply(df_pred, .(time, Religion), summarize, y=mean(values))
df_pred2$y.sd <- ddply(df_pred, .(time, Religion), summarize, y=sd(values))[, 3]
df_pred2$y.mean.sd <- ddply(df_pred, .(time, Religion), summarize, y=sd(values)/100)[, 3]

## PLOT 2 

df$time2 <- df_pred2$time2 <- c(1:288,  1:287)

df$Religion <- factor(df$Religion)
df_pred2$Religion <- factor(df_pred2$Religion)

levels(df$Religion) <- c("Christian", "Muslim")
levels(df_pred2$Religion) <- c("Christian", "Muslim")


####1.1  fitted nbinom values: duration ####

p <- mod_nb1$fitted.values
df$y.nbinom <- p

library(zoo)

df$y.rollmean <- unlist(tapply(df$y, df$Religion, rollmean, 
                               k = 15, fill = "extend"))
k <- c()

for (i in levels(df$Religion)) {
  k <- c(k, ksmooth(df$time2[df$Religion == i], df$y[df$Religion == i],
               "normal", bandwidth = 20)$y)
}

df$y.kernel <- k

df$facet_label <- ""

m <- tapply(df$y, df$Religion, mean)
means <- data.frame(mean = m, Religion = factor(levels(df$Religion)))


#######################################
#### 1.7 greyscale loess: duration ####
#######################################

p3_g <- ggplot(df, aes(x=time2, y=y)) + geom_point() + 
  stat_smooth(method = "loess", )
p3_g <- p3_g + facet_wrap(~Religion) + xlab("Rank") + ylab("Support for extremist statement from 1 (Totally disagree) - 7 (Totally agree)")
p3_g <- p3_g + theme(legend.position = "none") + theme_bw()
p3_g <- p3_g + ylim(c(-.1, 7.1)) + geom_hline(data = means, aes(yintercept = mean), linetype = "dashed")
p3_g

# combined plot: histogram
df$facet_label <- ""
p_hist_g <- ggplot(data = df, aes(x = y, fill = Religion)) + geom_bar(width=0.8, position = position_dodge(width=0.2))
p_hist_g <- p_hist_g + theme(legend.position = "none") + xlab("") + ylab("Count")
p_hist_g <- p_hist_g + coord_flip() + scale_fill_grey()
p_hist_g <- p_hist_g + facet_wrap(~facet_label)
p_hist_g <- p_hist_g + xlim(c(0.5, 7.5)) + theme_bw()
p_hist_g

# combined plot: empty:

empty_l <- ggplot()+geom_histogram(data = df, aes(x=time2, fill = Religion)) +
  theme(                              
    plot.background = element_blank(), 
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(), 
    panel.border = element_blank(), 
    panel.background = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank()
  ) + scale_fill_grey() + ylim(c(-5,0)) + theme(legend.position = "left") + theme(legend.title = element_blank())
empty_l <- empty_l + guides(col = guide_legend(reverse = TRUE))


# combined plot: all 

grid.arrange(p3_g, p_hist_g, 
             ncol=2, nrow=1, widths=c(2, 2), heights=c(1))

pdf("radicalization.pdf", width = 12, height = 6)
grid.arrange(p3_g, p_hist_g, 
             ncol=2, nrow=1, widths=c(2, 2), heights=c(1))
dev.off() 




